Minimizing Final Time (Jennings Problem)
Solve an optimal control problem with a minimal final time. Set up and solve the Jennings optimal control benchmark problem.
Problem Statement and Model
When solving differential equations over a variable time interval $[0,t_f]$, we can apply a time-scaling transformation to normalize the interval to$[0,1]$. This is achieved by introducing a final time parameter $t_f$. The Jennings optimal control problem divides derivatives by $t_f$. In practice, $t_f$ appears on the right hand side to avoid any divisions by 0.
\[\begin{gathered} \frac{\frac{dx}{dt}}{t_f} = f(x,u) \\ \frac{dx}{dt} = t_f f(x,u) \\ \end{gathered}\]
Our specific problem is defined as the following:
\[\begin{aligned} &&\min_{u(t),t_f} t_f \\ &&\text{s.t.} &&& \frac{dx_1}{dt}= t_f u, && t \in [0,1] \\ &&&&&\frac{dx_2}{dt} = t_f \cos(x_1(t)), && t \in [0,1] \\ &&&&&\frac{dx_3}{dt} = t_f \sin(x_1(t)), && t \in [0,1] \\ &&&&&x(0) = [\pi/2, 4, 0] \\ &&&&&x_2(t_f) = 0 \\ &&&&&x_3(t_f) = 0 \\ &&&&&-2 \leq u(t) \leq 2 \end{aligned}\]
Model Definition
First we must import $InfiniteOpt$ and other packages.
using InfiniteOpt, Ipopt, Plots;Next we specify an array of initial conditions.
x0 = [π/2, 4, 0]; # x(0) for x1,x2,x3We initialize the infinite model and opt to use the Ipopt solver
m = InfiniteModel(Ipopt.Optimizer);Recall t is specified as $\ t \in [0,1]$:
@infinite_parameter(m, t in [0,1],num_supports= 100)tNow let's specify descision variables. Notice that $t_f$ is not a function of time and is a singular value.
@variable(m, x[1:3], Infinite(t))
@variable(m, -2 <= u <= 2, Infinite(t))
@variable(m, 0.1 <= tf);Specifying the objective to minimize final time:
@objective(m, Min, tf);Define the ODEs which serve as our system model.
@constraint(m, ∂(x[1],t) == tf*u)
@constraint(m, ∂(x[2],t) == tf*cos(x[1]))
@constraint(m, ∂(x[3],t) == tf*sin(x[1]));Set our inital and final conditions.
@constraint(m, [i in 1:3], x[i](0) == x0[i])
@constraint(m, x[2](1) <=0)
@constraint(m, x[3](1) <= 1e-1);Problem Solution
Optimize the model:
optimize!(m)This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 1794
Number of nonzeros in inequality constraint Jacobian.: 2
Number of nonzeros in Lagrangian Hessian.............: 700
Total number of variables............................: 701
variables with only lower bounds: 1
variables with lower and upper bounds: 100
variables with only upper bounds: 0
Total number of equality constraints.................: 600
Total number of inequality constraints...............: 2
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 2
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0999999e-01 4.00e+00 7.38e-03 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.0009999e-01 3.98e+00 2.15e+02 -1.0 4.00e+00 - 1.00e+00 4.78e-03h 1
2 1.0000218e-01 3.98e+00 1.31e+05 -1.0 4.16e+01 - 4.77e-02 7.77e-05h 1
3 1.0000001e-01 3.98e+00 1.32e+05 -1.0 1.41e+05 - 1.41e-05 1.39e-05h 1
4 1.0000001e-01 3.98e+00 1.07e+05 -1.0 2.60e+04 - 9.52e-05 5.36e-05h 1
5 1.0000001e-01 3.98e+00 2.07e+05 -1.0 1.43e+04 - 1.90e-04 2.73e-05h 1
6 1.0000000e-01 3.98e+00 4.50e+05 -1.0 1.27e+05 - 2.01e-05 8.64e-06h 1
7r 1.0000000e-01 3.98e+00 1.00e+03 0.6 0.00e+00 - 0.00e+00 4.78e-07R 2
8r 1.0014680e-01 4.05e+00 1.92e+03 0.6 6.02e+00 - 6.22e-01 3.69e-02f 1
9r 1.0062971e-01 3.58e+00 1.41e+03 0.6 5.18e+00 - 1.00e+00 2.61e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10 1.0064745e-01 3.58e+00 3.90e+03 -1.0 3.99e+00 - 1.00e+00 2.60e-04h 1
11 1.0105387e-01 3.57e+00 6.14e+03 -1.0 8.03e+02 - 7.79e-04 4.93e-04h 1
12 1.0130844e-01 3.57e+00 2.30e+04 -1.0 3.32e+03 - 1.74e-04 4.12e-05h 1
13 1.0131275e-01 3.57e+00 2.13e+05 -1.0 1.09e+04 - 2.61e-04 3.20e-05h 1
14r 1.0131275e-01 3.57e+00 1.00e+03 0.6 0.00e+00 - 0.00e+00 4.88e-07R 5
15r 1.8100817e-01 3.50e+00 1.88e+03 0.6 4.36e+00 - 7.49e-01 5.12e-02f 1
16r 1.8669905e-01 2.33e+00 2.00e+02 0.6 2.29e+00 - 9.28e-01 7.89e-01f 1
17 1.8722777e-01 2.33e+00 1.68e+02 -1.0 3.81e+00 - 1.34e-01 8.11e-04h 1
18 1.8921350e-01 2.33e+00 2.25e+04 -1.0 1.06e+02 - 3.22e-02 2.37e-04h 1
19 1.9127731e-01 2.32e+00 3.87e+04 -1.0 4.33e+02 - 5.26e-03 3.02e-03h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20 1.9146732e-01 2.32e+00 1.58e+05 -1.0 4.59e+03 - 6.15e-04 1.52e-04h 1
21 1.9156826e-01 2.32e+00 9.58e+05 -1.0 8.57e+03 - 3.77e-04 7.65e-05h 1
22r 1.9156826e-01 2.32e+00 1.00e+03 0.4 0.00e+00 - 0.00e+00 3.51e-07R 5
23r 2.1757015e-01 1.96e+00 7.32e+02 0.4 1.33e+00 - 9.92e-01 2.67e-01f 1
24 2.1787418e-01 1.96e+00 3.86e+02 -1.0 3.24e+00 - 7.15e-02 1.89e-04h 1
25 2.1823093e-01 1.96e+00 9.27e+05 -1.0 2.86e+01 - 1.05e-01 4.33e-05h 1
26 2.1826491e-01 1.96e+00 4.64e+06 -1.0 2.91e+04 - 1.14e-04 2.25e-05h 1
27r 2.1826491e-01 1.96e+00 1.00e+03 0.3 0.00e+00 - 0.00e+00 1.55e-07R 2
28r 2.6610697e-01 1.14e+00 4.55e+02 0.3 1.52e+00 - 9.91e-01 5.45e-01f 1
29r 3.2355056e-01 1.72e-01 3.40e+02 0.3 4.55e+00 - 5.90e-01 2.47e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r 5.9791076e-01 2.47e-01 8.36e+01 0.3 1.08e+00 0.0 1.00e+00 8.10e-01f 1
31r 1.5393156e+00 1.30e+00 3.41e+02 -0.4 4.27e+00 0.4 1.78e-01 4.06e-01f 1
32r 1.7522485e+00 3.94e-01 1.40e+02 -0.4 1.26e+00 0.9 8.91e-01 6.62e-01f 1
33r 3.0683175e+00 1.80e+00 2.06e+02 -0.4 4.88e+00 0.4 2.65e-01 5.12e-01f 1
34r 4.0829865e+00 1.84e+00 2.89e+02 -0.4 1.66e+01 -0.1 2.27e-02 1.02e-01f 1
35r 4.2796802e+00 1.60e+00 3.08e+02 -0.4 2.03e+00 0.3 9.01e-01 1.26e-01f 1
36r 4.3283517e+00 1.51e+00 3.80e+02 -0.4 6.28e+00 - 1.00e+00 5.78e-02f 1
37r 4.3814442e+00 1.08e+00 4.96e+02 -0.4 5.56e-01 - 1.00e+00 2.82e-01f 1
38r 4.4392851e+00 1.45e-02 8.75e+00 -0.4 4.17e-01 - 1.00e+00 1.00e+00f 1
39 4.6154499e+00 1.30e-02 2.29e-02 -1.0 2.45e-01 - 9.73e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 4.4905204e+00 1.16e-02 2.94e-03 -1.7 3.27e-01 - 1.00e+00 1.00e+00h 1
41 4.3393561e+00 6.04e-02 1.38e-02 -3.8 9.67e-01 - 9.08e-01 9.34e-01h 1
42 4.2946439e+00 6.77e-02 8.52e-03 -3.8 1.80e+00 - 8.70e-01 1.00e+00h 1
43 4.2898253e+00 1.86e-02 1.26e-03 -3.8 2.46e+00 - 1.00e+00 1.00e+00h 1
44 4.2875206e+00 4.66e-03 1.80e-04 -3.8 2.82e+00 - 1.00e+00 1.00e+00h 1
45 4.2874424e+00 1.23e-04 3.37e-06 -3.8 6.32e-01 - 1.00e+00 1.00e+00h 1
46 4.2847486e+00 1.90e-03 3.96e-04 -5.7 1.13e+00 - 8.52e-01 1.00e+00h 1
47 4.2846175e+00 4.56e-04 2.07e-05 -5.7 1.34e+00 - 9.72e-01 1.00e+00h 1
48 4.2846039e+00 7.22e-05 5.56e-07 -5.7 6.40e-01 - 1.00e+00 1.00e+00h 1
49 4.2846034e+00 8.59e-07 6.40e-09 -5.7 5.87e-02 - 1.00e+00 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 4.2845650e+00 6.11e-06 1.18e-06 -8.6 1.22e-01 - 9.77e-01 1.00e+00h 1
51 4.2845649e+00 1.55e-08 1.17e-10 -8.6 1.04e-02 - 1.00e+00 1.00e+00h 1
52 4.2845648e+00 7.53e-12 5.91e-14 -9.0 1.52e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 52
(scaled) (unscaled)
Objective...............: 4.2845648348476271e+00 4.2845648348476271e+00
Dual infeasibility......: 5.9111156714079624e-14 5.9111156714079624e-14
Constraint violation....: 7.5290884637979616e-12 7.5290884637979616e-12
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 9.0923652519966649e-10 9.0923652519966649e-10
Overall NLP error.......: 9.0923652519966649e-10 9.0923652519966649e-10
Number of objective function evaluations = 67
Number of objective gradient evaluations = 41
Number of equality constraint evaluations = 67
Number of inequality constraint evaluations = 67
Number of equality constraint Jacobian evaluations = 57
Number of inequality constraint Jacobian evaluations = 57
Number of Lagrangian Hessian evaluations = 52
Total seconds in IPOPT = 0.191
EXIT: Optimal Solution Found.
Extract the results. Notice that we multiply by $t_f$ to scale our time.
ts = value(t)*value(tf)
u_opt = value(u)
x1_opt = value(x[1])
x2_opt = value(x[2])
x3_opt = value(x[3]);Plot the results
plot(ts, u_opt, label = "u(t)", linecolor = :black, linestyle = :dash)
plot!(ts, x1_opt, linecolor = :blue, linealpha = 0.4, label = "x1")
plot!(ts, x2_opt, linecolor = :green, linealpha = 0.4, label = "x2")
plot!(ts, x3_opt, linecolor = :red, linealpha = 0.4, label = "x3");Maintenance Tests
These are here to ensure this example stays up to date.
using Test
@test termination_status(m) == MOI.LOCALLY_SOLVED
@test has_values(m)
@test u_opt isa Vector{<:Real}
@test x1_opt isa Vector{<:Real}Test PassedThis page was generated using Literate.jl.